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Abstract. We have written a light curve synthesis code 
that makes direct use of model atmosphere specific inten- 
sities, in particular the NextGen model atmosphere grid 
for co ol gia nts (T e g < 6800 K and log (5) < 3.5, Hauschildt 



et al. 1999 ). We point out that these models (computed 



using spherical geometry) predict a limb darkening be- 
haviour that deviates significantly from a simple linear or 
two-parameter law (there is less intensity at the limb of the 
star) . The presence of a significantly nonlinear limb dark- 
ening law has two main consequences. First, the ellipsoidal 
light curve computed for a tidally distorted giant using the 
NextGen intensities is in general different from the light 
curve computed using the same geometry but with the 
black body approximation and a one- or two-parameter 
limb darkening law. In most cases the light curves com- 
puted with the NextGen intensities have deeper minima 
than their black body counterparts. Thus the light curve 
solutions for binaries with a giant component obtained 
with models with near linear limb darkening (either black 
body or plane-parallel model atmosphere intensities) are 
biased. Observations over a wide wavelength range (i.e. 
both the optical and infrared) are particularly useful in 
discriminating between models with nearly linear limb 
darkening and the NextGen models. Second, we show 
that rotational broadening kernels for Roche lobe filling 
(or nearly filling) giants can be significantly different from 
analytic kernels due to a combination of the nonspherical 
shape of the star and the radical departure from a simple 
limb darkening law. As a result, geometrical information 
inferred from V m t sin i measurements of cool giants in bi- 
nary systems are likewise biased. 
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1. Introduction 

The study of close binary stars is of interest for several 
reasons. For example, the understanding the structure and 
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evolution of stars is a basic goal of stellar astronomy, and 
is required in most other branches of astronomy. Critical 
tests of evolution theory (i.e. predicting the radius and 
luminosity of a star as a function of its mass and age) 
for stars other than the Sun are practical only for a small 
set of eclipsing bina ry stars (see , e.g. , Pols et al. 1997 ; 
Schroder et al. 1997 ; Lacy et al. |2000|) . In addition, our 
knowledge of the masses of stellar black holes (and many 
neutron stars and white dwarfs as well) depends on our 
ability to interpret ellipsoidal and usually non-eclipsing 
light curves (e.g. Avni fc Bahcall 1975; Avni 1978 



Clintock & Remillard p.990t Haswell et al. 
et al. 



1994; Sanwal et al. 1996; Orosz & Bailyn 



Mc- 
1993: Shahbaz 
Fi- 



1997) 



nally, considerable observational effort has been put forth 
recently to use detached eclipsing binaries as extragalac- 
tic distance indicators (Mochejska et al. 1998; Ribas et al. 
Since close binary stars are of such importance, it 



is crucial that we have the ability to construct accurate 
synthetic light curves for a variety of close binaries. 



The light curve expected from a particular close bi- 
nary depends on the system geometry (i.e. the figures of 
the stars, their relative sizes, separation, viewing angle, 
etc.) and on the radiative properties of the stars, which are 
set mainly by their effective temperatures, surface gravi- 
ties and chemical compositions. The equations describing 
the basic system geometry for a close binary are reason- 
ably simple and have been known for a long time (e.g. 
see the text by Kopal 1959). In practice, however, the di- 
rect computation of light curves requires a computer, and 
codes to compute light curves have been in use since the 



late 1960s (e.g. Lucy |967] Hill & Hutchings 1970} Wil- 
son & Devinney 1971; Mochnacki & Doughty |1972|; Avni 
& Bahcall 1975 ; see also the review by Wilson 1994 ). On 
the other hand, the detailed computation of stellar atmo- 
sphere models which describe the specific intensity of radi- 
ation emitted by the stellar surfaces is quite involved. As 
a result, approximations are frequently used in the com- 
putations of light curves. The Planck function is usually 
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used to compute the normal monochromatic intensity of 
a surface element with an effective temperature T e ff. 



Iq oc [exp(ft,c/fcAr c ff) — 1] 



(1) 



where A is the effective wavelength of the observation and 
h, c, and k are the usual physical constants. After Iq is 
found, the specific intensities for other emergent angles 
are computed using a simple parameterisation (the "limb 
darkening law" ) : 



/(/i)=/ [l-x(l-/i)]. 



(2) 



The coefficient x depends on the temperature, gravity, and 
chemical composition of the star being modelled. There 
are also many two-parameter limb darkening laws includ- 
ing a "quadratic" law: 

/( M ) = /o[i-z(i-M)-y(i-/i) 2 L (3) 

the "square root" law (Di'az-Cordoves & Gimenez f 992| ): 

I(fi)=I [l-x(l- t i)-y(l-^JI)], (4) 

or the "logarithmic" law (Klinglesmith & Sobieski f 97C| ): 

= -M 1 ~ x{\ - fi) - y/iln/i]. (5) 

There are many tabulations of the coefficients for these 
various laws for a wide variety of temperatures, gravities, 
and chemical compositions (e.g. Al-Naimiy 1978; Wade & 
Rucinski |1985|; Van Hamme |1993j; Claret |1998|). 



We discuss in this paper our technique for computing 
light curves using model atmosphere specific intensities di- 
rectly, thereby eliminating the need for the black body ap- 
proximation and the one- or two-parameter limb darken- 
ing laws. Although our emphasis here is on stars with low 
effective temperatures and surface gravities (T e g < 6800 K 
and log g < 3.5), the technique is quite general. In the next 
section we discuss some previous work in this area and give 
details of our method. We then discuss some of the basic 
results and their implications, and give a summary of this 
work. We have also included an appendix to the paper 
which gives some of the details of our light curve code not 
directly related to the stellar atmospheres. 

2. Computational Technique 

2.1. Previous Work 

The use of model atmosphere computations inside a light 
curve synthesis co de is b y no means new. The widely used 
Wilson-Devinney ( 1971 ; hereafter W-D) code has routine 
which applies a correction to the normal (black body) in- 
tensity: 



rcorr 
J 



I r(T eS ,logg, A), 



(6) 



where r(T e ff, log g, A) is the ratio of the filter-integrated 
stellar atmosphere model characterised by T e g and log g 
to the filter-integrated blackbody intensity. The specific 



intensities for other angles are then computed from 7™"' 
using the limb darkening law. The correction routine sup- 
plied with the W-D code is based on the Carbon & Gin- 



gerich (1969) models. R. E. Wilson (priv. comm.) informs 
us he is in the process of updating the correction rou- 



tine. Milone et al. (1992) have independently written a 
correction routine for the W-D code based on the Kurucz 
( |1979[ ) models. Linnell & Hubeny Q1994 |199€| ) have writ- 
ten a series of codes to compute synthetic spectra and 
light curves for binary stars, including ones with disks. 
They use Hubeny's ge neral spectrum synthesis code SYN- 
SPEC (Hubeny et al. 1994|) t o generate the model spectra, 
using as input Kurucz ( 1979] ) atmosphere models for cool 
stars (T efi < 10000 K) and T LUSTY (Hubeny |l988l ) and 
TLUSDISK (Hubeny |l99l[ ) models for hotter stars and 
disks, respectively. They have applied their model with 
some success to (3 L yrae (L innell et al. 1998a|) an d MR 
Cygni (Linnell et a l. |l998b|) . Tjemkes et al. ( |l986|) have 
used Kurucz ( 1979] ) models in their light curve synthe- 
sis code by tabulating filter-integrated normal intensities 
for a grid of effective temperatures and surface gravities. 
Specific intensities for other emergent angles are computed 
from a limb darkening law tabulated from Kurucz models. 
This code has been applied successfully to X-ray binaries 
such as LMC X-4 (Heemskerk & van Paradijs 1989) and 
GRO J1655-40 (van der Hooft et al. |1998| ) and to th e 
millisecond pul sar P SR 1957+20 (Callanan et al. |1995[) 
Shahbaz et al. (1994) mention of the use of Kurucz (1979) 
fluxes in a code used at Oxford, but this paper does not 
specifically describe how the model atmosphere fluxes are 
incorporated i nto th e light curve synthesis code. Similarly, 
Sanwa l et al. ( 1996 ) mention the use of Bell & Gustafsson 
( |l989|) model atmosphere fluxes in the light-curve synthe- 



1986 



Haswell et al. 1993) without giving 



sis code devel oped at the Universi ty of Texas at Austin 
(Zhang et al. 
specific details 



2.2. Current Work 

We use a technique suggested to us by Marten van Kerk- 
wijk to incorporate model atmosphere intensities into our 
light curve synthesis code. In general, a detailed model at- 
mosphere computation yields the specific intensity /(A, /i) 
at a given wavelength A and emergent angle fi = cos 9, 
where fj, = 1 at the centre of the apparent stellar disk and 
/j = at the limb. The total disk-integrated intensity /(A) 
observed at the wavelength A is then 



/(A) = / /(A.mWm- 
Jo 



(7) 



Normally, the light curves of a binary star are observed in 
a broad bandpass. The observed intensity Tfilt hi a given 
bandpass is 



FILT — 
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I(X, p,)/jLdfj, 



W FIhT (X)dX, 



(8) 



where Wfilt is the wavelength-dependent transmission of 
the filter bandpass in question. We can reverse the order 
of the integrals in Eq. (@) to give 



h 



FILT 



I(X, f i)W F1LT (X)dX 

7filt(^)MM' (9) 
The quantity in square brackets in Eq. (|^), namely 

r+oo 

/filt(m)= / I(X,ti)W FlLT (X)dX (10) 



is independent of any geometry (for non-irradiated atmo- 
spheres), and as such can be computed in advance of the 
light curve synthesis computations. 

In our current implementation we compute for each 
model characterised by a temperature T g and gravity 
log g a table of eight filter- integrated intensities for each 
specific angle /i. The eight filters currently are the Johnson 
U BVRIJHK system where we have used the optical filter 



response curves given in Bessell ( L990) and the infrared fil- 
ter response curves given in Bessell & Brett ( 1988| ). Using 



all of the models we then generate a table of the form 



Tl 


gl 








Nmu 










mul 


1(1,1) 


1(1,2) . 


. 1(1,7) 


1(1,8) 


mu2 


1(2,1) 


1(2,2) . 


. 1(2,7) 


1(2,8) 


muN 


I(N,1) 


I(N,2) . 


. I(N,7) 


I(N,8) 


Tl 


g2 








Nmu 










mul 


1(1,1) 


1(1,2) . 


. 1(1,7) 


1(1,8) 


mu2 


1(2,1) 


1(2,2) . 


. 1(2,7) 


1(2,8) 


muN 


I(N,1) 


I(N,2) . 


. I(N,7) 


I(N,8) 


T2 











The first two numbers are the temperature and gravity of 
the first model. The next line gives the number of specific 
intensities that follow. The next iV M lines give the value 
of n followed by the eight filter-integrated intensities. The 
format is then repeated for the additional models. The 
table is sorted in order of increasing temperature, and for 
each temperature, the table is sorted in order of increasing 
gravity. 

During the course of computing a synthetic light curve 
the specific intensity of each surface element must be spec- 
ified. These values of I(T U i P , log <?i n p , Minp) are interpolated 
from the table using a simple linear interpolation proce- 
dure. First, we locate the 4 nearest models (T up , logg up ), 
(T up ,loggi ow ), (Tiow^oggup), and (T iow , log3i ow ), where 



T eff =4000, 
• V 
» I 


i 1 i 1 i 1 i 

log(g) = 1.0 


J- t 

■ m-J-\ — 




; T eff =5000, 

• u 


-H 1 1 1 1 1 1 ; 

log(g) = 0.0 " 


'. * K 




-/ 




"* . 

i i ilnn» 


wd&r 1,1,1, 



Fig. 1. Top: the filter-integrated specific intensities 
-^filt(m) f° r Johnson V and I for a model with T c g = 
4000 K and log<? = 1. The dashed-dotted lines show the 
limb darkening behaviour predicted by the "square root" 
law (Eq. (|])), where we have used coefficients (computed 
from Kurucz models) given in Van Hamme (1993). Bot- 
tom: the filter-integrated specific intensities /filt(a*) for 
Johnson U and K for a model with T e g = 5000 K and 
\ogg = 0, and the corresponding square root limb darken- 
ing laws. 



?low < Tinp < T nv and loggi ow < logg inp < log5 up . Next, 
within each of the 4 nearest models we find the filter- 
integrated intensities appropriate for fJ,i np , using linear in- 
terpolation. Finally, the desired values of 
7(Ti np , log(/i n p, /Xinp) are found by linear interpolation first 
in the log g direction and then the T c g direction. Returning 
the filter-integrated intensities for eight different niters at 
once is advantageous when modelling observations in sev- 
eral bandpasses since one has to search the intensity table 
only once per surface clement per phase. 

As we noted above, the technique we just outlined is 
quite general. We now turn our attention to the recent 
grid of NextGen models for cool giants computed by 
PHH using the parallelised version of his general purpose 



code PHOENIX (Hauschildt et al. p_997| ). The computa- 
tional techniques and the input physics used to compute 
the NextGen models for c ool gi ants have been discussed 
elsewhere (Hauschildt et al. 1999 ) and will not be repeated 
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'g(g) (cgs) 



Fig. 2. Contour plots of the quantity /i cut for the 
NextGen models in the B (top), I (middle), and K (bot- 
tom) bands. There are no static models for the highest 
temperatures and lowest gravities, hence there are gaps in 
the upper left parts of the contour maps. 




12 3 



log(g) (cgs) 

Fig. 3. Contour plots of the "lost light," expressed as a 
per cent, of the NextGen models in the B (top), I (mid- 
dle), and K (bottom) bands. There is a large area of the 
parameter space for the B filter (shown as the large gap 
in the contours) where there is no lost light as defined in 
the text. 



here. These spherical models require an input mass to 
compute the sphericity factor, and we use M = 5M Q 
for all of the models. We did examine a few models with 
masses of 2.5 Mq and 7.5 Mq and found no discernible dif- 
ferences compared to the M — 5 Mq models. For the mo- 
ment we restrict ourselves to the models with solar metal- 
licity. For each model the specific intensities are computed 
for 64 different angles over a wavelength range of 3000- 
24,998 A in 2 A steps. The distribution of the emergent 
angles fj, is chosen by the PHOENIX code based on the 
structure of the atmosphere under consideration. Hence, 
in general, the 64 values of /i are different from model to 
model. 

Fig. |l| shows the quantity /filt(m) f° r two different 
NextGen models and filter combinations. The top panel 
shows the integrated intensities in the Johnson V and / fil- 
ters for the model with T c s — 4000 K and logg = 1.0 and 
the bottom panel shows the integrated intensities in the 
Johnson U and K filters for the model with T e g = 5000 K 
and logg = 0.0. For comparison, we also show the square 
root limb darkening laws (Eq. (^)) computed usin g coe f- 
ficients taken from the tabulation of Van Hamme ( 1993[ ). 
These four curves essentially show the limb darkening be- 



haviour predicted by the Kurucz models to within a few 
per cent, although is important to note that the intensi- 
ties for /i = 1.0 are different between the Kurucz models 
and the NextGen models (see Hauschildt et al. (1999) 
for further discussion of this point). This figure is quite 
striking. It shows the radically different limb darkening be- 
haviour predicted by the spherical NextGen models and 
the plane-parallel Kurucz models. The NextGen model 
predicts a sharp decrease in the intensity at relatively large 
fj, values (as large as ss 0.4), whereas the limb darken- 
ing parameterisations predict substantial intensity all the 
way to fi = 0. This sudden decrease in the intensity is 
a consequence of the spherical geometry: for sufficiently 
low gravities (log<? 3.5 for most temperatures) there is 
much less material near the limbs and hence much less ra- 
diation. The radiation that would have come out near the 
limb in the redder bandpasses comes out at much shorter 
wavelengths (see Hauschildt et al. 1999 ). 

The value of \i where the sudden fall-off in the intensi- 
ties occurs depends on the effective temperature and grav- 
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ity of the stellar model, and on the bandpass. We define 
the slope of the intensities at to be 



s(i) 



J(i + 1) - J(i - 1) 
M(* + i) - /*(* - 1) 



(11) 



The location /i cu t of the edge is defined to be the point 
where s(i) is the largest. Fig. [2] shows contour plots of /i cu t 
for three different bands (Johnson B, I, and K). In gen- 
eral, for a fixed effective temperature, /i cu t gets larger as 
the gravity of the model decreases. Also, for a fixed log g, 
the hotter models tend to have a slightly larger values of 
// cu t- There is some irregular behaviour seen in the con- 
tours of ii cu t at temperatures between about 3000 K and 
4000 K. The reason for this seems to be that for models 
in this temperature range, the fall-off in the intensity is 
generally less sudden, especially in the bluer filters. Thus 
/i cu t is less clearly defined in these cases. 

It is quite obvious that the quantity /filt (Eq. (0)) 
evaluated for the NextGen intensities will be less than 
the corresponding quantity evaluated for a one- or two- 
parameter limb darkening law. For each model and filter, 
we computed the "lost light" in the following way. We first 
"linearised" our table of NextGen intensities: We fitted 
a line to the intensities for the ten specific angles nearest 
to /i = 1.0. Then the specific intensities for the remain- 
ing 54 angles were replaced by the intensities determined 
from the extrapolation of the fitted line. In most cases, 
the extrapolated intensities at the limb were above zero. 
There were a few cases such as the T e ff = 5000, log g = 0.0 
U -band model shown in Fig. [I] where the extrapolated in- 
tensities at the limb would have been negative. In such 
cases we replaced all negative intensities with zero. Eq. 
@ was evaluated for the linearised table and the regu- 
lar table, and we define the lost light as D(T e s, log g) = 

100 X (/FILT.linearPrfF^Oga) - iFILT,rcgular (T c ff , log g) ) / 

^FiLT,iincar(Tcff,logsO. Fig. j| shows contour plots of 
D(T c ff , log g) for the B (top), I (middle), and K (bottom) 
filters. By our definition the lost light in many of the B- 
band models is zero. Otherwise, for the B band, the lost 
light is typically a few per cent and about 3 or 4 per cent 
for the hottest models. On the other hand, the lost light in 
the / and K bands can be quite substantial, up to 17 per 
cent for the models with lowest gravities and the hottest 
temperatures. 

We can easily make an image of the disk a model star as 
it would appear in the sky. The intensity maps displayed in 
Fig. H compare the NextGen intensities for T cS = 5000 K 
and log g = 0.0 to monochromatic black body intensities 
(limb darkened using the square root law (Eq. (Q)) for 
three different filters (Johnson B,V, and I). The star with 
the NextGen intensities appears smaller on the sky. The 
difference in the limb darkening closer to the disk centre 
is also apparent from the figure. 

Finally, for comparison purposes, we generated an in- 
tensity table from Kurucz models. We retrieved the file 
"ip00k2.pckl9" from Kurucz's public World Wide Web 



pagesfj. This file contains the specific intensities for 17 dif- 
ferent angles for a grid of solar abundance models. The 
coolest models have T e g = 3500 K. All of these models 
are plane-parallel and LTE. 

3. General Results 

We have written a new light curve synthesis code which 
is based on the work of Avni (Avni & Bahcall 1975; 
Avni 1978). This code, named ELC, can model any semi- 



detached or detached binary with a circular orbit (we will 
generalise the code to include overcontact binaries and 
binaries with eccentric orbits within the next year). The 
second star can be surrounded by an accretion disk. We 
will defer giving specific details of the code to the Ap- 
pendix and now turn our attention to the results related 
to the use of the NextGen intensities. 

3.1. Light curve amplitude 

For various test cases we computed light curves in four dif- 
ferent ways: using our ELC code in the black body mode, 
using the W-D code in the black body mode (we have 
added a phase shift of 0.5 to the W-D models, see the Ap- 
pendix) , using our ELC code with the Kurucz intensity ta- 
ble, and using our ELC code with the NextGen intensity 
table. Fig. || shows four different normalised light curves 
(in V, I, and K) for a single Roche lobe filling star in 
synchronous rotation with T c s = 4000 K and logg = 1.0. 
The mass ratio is 5, meaning the unseen second star is 
five times more massive than the cool giant. The incli- 
nation is 75 degrees. For the black body models we used 
the square ro ot lim b darkening law with coefficients from 
Van Hamme ( 1993 ). The gravity darkening exponent was 
f3 = 0.08, appropriate for a star with a convective en- 
velope. (We note that the black body light curves were 
computed using a single limb darkening law for the en- 
tire star, whereas the light curves computed using model 
atmosphere intensities (NextGen or Kurucz) will have 
location-specific limb darkening "built in".) In all three 
cases, the two black body curves are nearly identical (the 
largest deviations between the two are on the order of 2 
millimags). Typically the NextGen curves have deeper 
minima than the black body curves. Depending on the fil- 
ter, the NextGen light curves are different from the the 
light curves computed with Kurucz intensities, although 
for the example shown in Fig. ^ the two /-band curves 
have only minor differences. 

In general, the light curves computed using the Next- 
Gen intensity table are different from light curves com- 
puted using the same geometry but with black body in- 
tensities. This comes as no surprise, given the potentially 
large differences in the limb darkening between the Next- 
Gen models and the one- or two-parameter limb darken- 
ing laws. It is usually the case (especially in the /, J, 

1 http://cfaku5.harvard.edu/ 
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T p „=5000, log(g)=0 B 



blackbody 
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0.4 
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Fig. 4. Intensity maps for a spherical star with T c g = 5000 K and log g = 0.0. The top images show intensities computed 
from the NextGen models in the Johnson B, V, and / bandpasses, while the bottom images show monochromatic 
black body intensities with the square root limb darkening law (Eq. (|])) at the effective wavelengths of the same three 
filters. In each case the intensities have been normalised to 1.0 at the disk centre. 



H, and K filters) that the NextGen light curves have 
larger amplitudes (in normalised intensity or in magni- 
tudes) than the corresponding black body light curves. 
Fig. shows why this is so. The top curve is the /-band 
light curve of the single star described in the Fig. |5| caption 
made with the linearised intensity table. Since the square 
root limb darkening law for this model (T e s = 4000, 
log<? = 1.0) and bandpass closely matches the actual 
NextGen intensities near disk centre (see Fig. [l]), the 
corresponding black body light curve is nearly identical 
to the "linearised" light curve. The middle curve in Fig. 
H is the light curve for the same star made with the reg- 
ular NextGen intensity table, plotted on the same in- 
tensity scale. The average light level in this light curve is 
about 5 per cent lower than the linearised one, which is 
roughly what we would expect from Fig. |[ The bottom 
curve is the difference between the two top light curves. 
This lower curve is basically the light lost near the limb, 
and it is essentially constant with phase. The two top 
curves have roughly the same amplitude in these intensi- 



ties units. That is, /„ 



is the same for both curves. 



However, the relative amplitude (I max — I m i n ) / I&ve for the 
regular NextGen light curve will be larger since its mean 
light level is lower. Hence, the amplitude of this particu- 
lar NextGen light curve (and most others) in normalised 
intensity or in magnitudes is larger than the amplitude of 
the corresponding black body curve. 



There are some situations such as the U model shown 
in the bottom of Fig. |l] where the square root limb darken- 
ing law does not match the slope of the NextGen inten- 
sities near fj, = 1. For cases like this, the black body light 
curve can have a larger amplitude than the NextGen 
light curve. 

Fig. H nicely shows the difference in the light curves 
between black body intensities and model atmosphere in- 
tensities (e.g. either NextGen or Kurucz). However, in 
the example shown, the NextGen light curves and the 
Kurucz light curves had only minor differences. To illus- 
trate how large the difference between a NextGen light 
curve and a Kurucz light curve can be we consider a binary 
similar to RZ Scuti. This semidetached binary consists of 
a B3Ib star which rotates near its centrifugal limit and 
a Roche lobe- filling A2? star (e.g. Olson & Etzel |1994| ). 
The mean gravities of the two stars are loggi « 3.2 and 
log 32 ~ 2.4, so sphericity effects should be important here. 
As a result of its rapid rotation near the critical limit (see 
Van Hamme & Wilson 199C for a discussion of "critical ro- 



tation"), the B-star is significantly flattened. Furthermore, 
the surface gravity near its equator is relatively low, giv- 
ing rise to a relatively large range in temperatures owing 
to the increased gravity darkening. We currently cannot 
model the real RZ Scuti binary using the NextGen inten- 
sities because the mean temperature of the B-star is well 
outside our current model grid. Indeed, there are parts 
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ELC, W-D black body 
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0.4 0.6 0.8 1 1.3 
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Fig. 5. Top: F-band light curves for a single Roche lobe 
filling star computed in four ways. The mean temperature 
of the star is T\ — 4000 K, its surface gravity is log g = 1.0, 
the binary mass ratio is Q = 5, and the inclination is i = 
75 degrees. The y-axis is an intensity scale. The two curves 
with the largest intensities at phase 0.5 are the ELC and 
W-D monochromatic curves computed using black body 
intensities (the two curves are nearly identical). The curve 
in the middle drawn with the dash-dotted line is the light 
curve generated by ELC using the Kurucz intensity table. 
The curve with the lowest intensity at phase 0.5 is the ELC 
curve computed using the NextGen intensities. Middle: 
Same as the top, but for the / filter. In this case the curves 
computed using the Kurucz and NextGen intensities are 
nearly the same. Bottom: Same as the top, but for the K 
filter. The ELC NextGen curve has the deepest minimum 
at phase 0.5. 



on the B-star where the temperature and gravity combi- 
nation fall outside the Kurucz grid, so we cannot model 
the real binary with the Kurucz grid either. For this ex- 
ample we therefore modified the temperatures of the two 
stars and "slowed down" the mass gaining primary slightly 
so that the temperature and gravity combination of each 
point on each star is contained in both the NextGen and 
Kurucz grids. The adopted model parameters are sum- 
marised in the caption of Fig. ^. Fig. [j] itself shows the 
difference between the NextGen light curve and the Ku- 
rucz light curve in magnitudes for three filters. There are 
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Fig. 6. The top curve is the /-band light curve for a single 
Roche lobe filling star with T e s = 4000 K and a mean 
gravity of log g = 1.0 constructed using the "linearised" 
intensity table. As in Fig. |5[ i = 75° and Q — 5. The I- 
band black body light curve shown in Fig. |5| is overplotted 
with the dash-dotted line (there is essentially no difference 
between the two curves). The middle curve is the light 
curve made using the regular NextGen intensity table, 
plotted on the same intensity scale. The bottom curve is 
the difference between the linearised curve and the regular 
curve. 



large differences between the two models, and the size of 
the difference depends on the filter bandpass. The differ- 
ence between if-band light curves is as large as 0.05 mag 
near phase 0, when the cooler star is eclipsed by the hot- 
ter star. One can also note from Fig. [7] that the U and V 
difference curves have means which are less than zero. In 
other words the binary is bluer in U — K and V — K when 
NextGen intensities are used compared to when Kurucz 
intensities are used. 

To quantify how much the difference in the light curve 
amplitudes between the NextGen models and the black 
body models might matter when fitting the light curves 
of a real binary star, we consider a model binary simi- 
lar to T Coronae Borealis. T CrB is an S-type symbiotic 
binary where an M4 giant orbits an optically faint hot 



companion with a period of about 227.6 days (Kraft 1958 
Kenyon 



Garcia 1986 



2000) 



There are no 
1992}. Never- 



Fekel et al 
eclipses observed in the UV (Selvelli et al 
thcless, the large amplitude of the ellipsoidal light curves 
suggests the M giant fills its Roche lobe and is viewed at 



a large inclination (Bailey 1975| ; Lines et al. 1988; Yudin 



Munari |1993Q . Since the density of a Roche lobe fill 



ing star is a function only of the orbital period to a good 



approximation (Faulkner et al. 1972; Eggleton 1983), the 



surface gravity of the M giant is only a weak function of 
its assumed mass. In the case of T CrB, \ogg w 0.7, so 
we expect the sphericity effects to be important her e. The 
mass ratio of the binary is not well known. Kraft ( 1958 ) 
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Fig. 7. The difference curves in magnitudes between the 
NextGen light curves and the Kurucz light curves for a 
binary with a geometry similar to RZ Scuti (i.e. a Roche 
lobe-filling "cool" star and a "hot" mass-gaining star ro- 
tating near its critical limit), but with altered tempera- 
tures. The adopted parameters are: i — 85.65°, Q — 0.21 6, 
P = 15.2 days, a = 61.89 R e , T x = 6000 K, T 2 = 5000 K, 
fx = 0.66045, h = 1.0, fii = 4, Q 2 = 1, ft = 0.25, 
ft = 0.08, and detailed reflection with albedos of 1.0 and 
0.5 for the hotter and cooler star, respectively (see the 
Appendix for a detailed discussion of the free parameters 
of the model). The resulting physical masses, radii, and 
mean gravities of the two stars are Mi = 11.3 Mq, R\ = 
14.4 R Q , log g 1 = 3.17 and M 2 = 2.45 M , R 2 = 16.2 R Q , 
\ogg 2 — 2.42, respectively. 



measured hydrogen emission line radial velocities on seven 
plates and found M g ; ant /M comp ss 1.4. However, this mass 
ratio implies a rotational velocity of V m t sin i « 23 km s _1 , 
which is much larger than the upper limit of ss 10 km s _1 
measured by Kenyon & Garcia (1986). We discuss below 
some of the potential problems associated with measuring 
the rotational velocity of a tidally distorted cool giant, 
and in view of this discussion the rotational velocity up- 



per limit of Kenyon & Garcia (1986) should be treated 
with caution. 

The J and K band light curves of T CrB collected be- 
tween 1987 August and 1995 June are stable and represen- 
tative of the true ellipsoidal modulation (Yudin & Munari 
1993j; Shahbaz et al. [19971). The ful1 amplitude of the J 



light curve is about 0.2 magnitudes. On the other hand, 
the V light curve shows additional sources of variability 
not associated with the underlying ellipsoidal modulation 
(e.g. Lines et al. 1988] ). Nevertheless, between about 1989 



and the beginning of 1996, the V light curve seemed to be 



reasonably stable. Hric et al. (1998) refer to this period 
as the "quiet stage". The full amplitude of the V light 
curve during this quiet stage is about 0.4 magnitudes. For 
this discussion we will assume this light curve represents 
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Fig. 8. The results of fitting a NextGen light curve using 
black body intensities. The input model parameters are 
single Roche lobe-filling star, T c g = 3560 K, Q = 1.6666, 
and i = 60°. Top panel: the O — C residuals (in magni- 
tudes) of a fit to the J-band light curve only (Q fixed, i 
free). The fitted inclination is i = 65.15°, and the ampli- 
tude of the 1^-band light curve is nearly 0.03 magnitudes 
too small. Middle panel: similar to top, but a fit to both 
V and J. The fitted inclination is nearly 11° too large, 
and the V residuals still have relatively large systematic 
deviations. Bottom panel: similar to the middle, but both 
the inclination and mass ratio were free parameters. The 
residuals are reasonably small, but the fitted parameters 
are quite different from the input parameters. 



the true ellipsoidal light curve. Both Shahbaz et al. (1997) 
and Belczynski & Mikolajewska ( 1998| ) had difficulty fit- 
ting the optical (V or /) light curves simultaneously with 
the infrared (J or K) light curves. If the amplitude of the 
J band light curve was matched, then the model V and / 
light curves (computed using black body intensities) had 
amplitudes that were about a factor of two too small. 

A thorough analysis of existing T CrB data is beyond 
the scope of the present paper and will be deferred to a 
future paper. For now it will suffice to discuss simulated T 
CrB-like light curves. For these light curves we will use the 
T CrB sy stem parameters adopted by Belczynski & Miko- 
lajewska ( 1998 ), namely M g i ant /M comp = 0.6, which in our 
notation is Q = 1/0.6 = 1.6666, i = 60°, T cS = 3560 K, 
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and a gravity darkening exponent of j3 = 0.08, appropri- 
ate for a star with a convective outer envelope. We also 
used the orbital period and if-velocity of the M giant 
given in Kenyon & Garcia ( |1986| ). Using the NextGen 
intensity table we computed light curves for the Johnson 

V and J filters. The model light curves were converted 
to magnitudes for compatibility with our optimising rou- 
tines. We used ELC in its black body mode to fit various 
combinations of the simulated light curves, and the re- 
sults are shown in Fig. ||. The top panel shows the results 
of fitting the J light curve only, using the inclination as 
the sole free parameter. The fitted inclination is 65.15°, 
and the O — C residuals for the J band are reasonably 
small 0.01 mag). However, if we take the predicted 

V light curve for this geometry and compare it to the 
simulated V light curve, we find that the amplitude of 
the black body V light curve is much smaller than the 
amplitude of the simulated V (NextGen) light curve. 
Thus we have reproduced the same problem that Shah- 
baz et al. ( |1997| ) and Belczyhski & Mikolajewska ( [1998D 



had when they attempted simultaneous optical and in- 
frared fits. Not surprisingly, if we attempt to fit both our 
simulated V and J light curves simultaneously using the 
inclination as the only free parameter, the fits to the J 
light curve get worse (middle panel). Finally, we fit both 
light curves simultaneously using both the inclination and 
mass ratio as free parameters. The O — C residuals for 
both filters are not excessively large and might be compa- 
rable to the observational errors in the real binary. How- 
ever, the fitted parameters are quite different from the 
input ones: Qm = 2.41 and igt = 67.2°, compared to in- 
put values of Qi np = 1.67 and ii np = 60.0°. Naturally, the 
derived component masses from the fit are significantly 
different from the input masses: M g i ant (fit) = 0.32 M Q 
and M comp (fit) = 0.78 M©, compared to input values of 
M giant (inp) = 0.71 M Q and M comp (inp) = 1.18 M Q , 



3.2. Rotational broadening kernels 

Measurements of the rotational velocities of the stars in 
close binaries provide powerful constraints on the light 
curve solution. Indeed, any light curve solution that spec- 
ifies the mass ratio, inclination, and the angular velocity 
ratios predicts specific values for the observed values of 
V m t sin i, provided of course that at least one radial veloc- 
ity curve is available. In cases where a star fills its critical 
lobe exactly and is in synchronous rotation (generally safe 
assumptions in cataclysmic variables and low mass X-ray 
binaries), a measurement of its rotational velocity con- 
strains the mass ratio of the binary (e.g. Wade & Horne 



1988 



If spectra with high resolution and high signal-to-noise 
are available, then one can use Fourier techniques to mea- 



sure the rotational velocity (e.g. Gray 1992 and cited ref- 
erences). When only spectra of lower quality are avail- 
able, then often one measures K-otsini by comparing the 

















, : 

■ . ■ 


~* " ^ XV 

■ x \ 


_ 








■ "\ 






r ■ 




• \ 




- 


f 






- 




ft 




















/ / 
























\ V 






' / ■ // 








J, 












//- 










// 




\\ 






,if 




Vv 






■If 




Vv. 





-1 



-0.5 

AA/X L 



0.5 



Fig. 9. Rotational broadening kernels G(A) for the star de- 
scribed in the caption to Fig. || for three different phases: 
0.0 (lower curves), 90° (middle curves, offset by 0.5), and 
180° (upper curves, offset by 1.0) (the giant is behind its 
invisible companion at phase 180°). The solid lines are 
the kernels for Roche geometry and NextGen intensities. 
The dash-dotted curves are for Roche geometry and black 
body intensities plus a linear limb darkening law with a 
coefficient of x = 0.6. The dotted curves are analytic ker- 
nels with a limb darkening coefficient of 0.6. 



spectrum of interest with a spectrum of a slowly rotat- 
ing template star (observed with similar instrumentation) 
that has been convolved with a broadening kernel G(A). 
Various trial values of the width of the broadening kernel 
are tried until the optimal match is found. Gray ( |1992| , 
pp. 370-374) gives a clear description of how to compute 
the broadening kernel G(A) analytically for the case where 
the intrinsic line profile H(\) has the same shape over the 
entire disk. Essentially, one can place an x,y coordinate 
system on the apparent disk of the star on the sky where 
the y-axis is parallel to the axis of rotation. The disk of 
the star can then be divided up into a number of strips 
parallel to the y-axis, each having its own Doppler shift 
according to its x-coordinate. The broadening kernel is 
evaluated at a particular Doppler shift by integrating the 
intensity over the appropriate strip. Normally one assumes 
a linear limb darkening law with a coefficient of x — 0.6 
when computing G(A). Finally, the broadened line pro- 
file is the convolution of the intrinsic line profile with the 
kernel: H Iot (X) = G(A) * H{\). 

The above discussion applies to spherical stars. A star 
that fills a substantial fraction of its Roche lobe departs 
markedly from spherical symmetry, and as such will have 



distorted line profiles (e.g. Kopal 1959: Marsh et al. 1994; 
Shahbaz 1998| ). Furthermore, the degree of the line profile 
distortion depends on the Roche geometry (i.e. mass ra- 



tio and inclination) and the orbital phase (Shahbaz 1998). 
The intensity maps that the ELC program can write (e.g. 
Fig. H) can be used to numerically evaluate G(A), thereby 
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Fig. 10. Rotational broadening kernels G(\) for UU Can- 
cri. The solid line is the kernel computed for a mass ratio 
of Q = 1.65, which corresponds to a "mean" rotational 
velocity of 27ri? eff /P = 20.7 km s _1 . The dash-dotted line 



is the analytic kernel computed by Eaton et al. ( 1991 ) 
corresponding to Q — 1.235 and 2ttR c h/P — 25.0 km s _1 



accounting for both the nonspherical shape of the star 
and the nonlinear limb darkening. Fig. ^| shows rotational 
broadening kernels G(A) for the star described in the cap- 
tion to Fig. H for phases 0.0, 90°, and 180°, when the giant 
is behind its invisible companion. The situation is compli- 
cated. Near phase 0.0, the projected surface of the giant 
in the plane of the sky is not too different from a circle 
(it is flattened slightly in the y-direction). Thus the ob- 
served profile will be nearly symmetric. However, the ker- 
nel, computed using the NextGen intensities, is narrower 
than the analytic kernel and the kernel computed using 
black body intensities and linear limb darkening since the 
star with the NextGen intensities appears smaller on the 
sky (Fig. ||). At phase 90° (quadrature), the star is elon- 
gated in the x-direction. As a result, the broadening kernel 
computed from the black body intensities is wider than 
the analytic kernel, and it is asymmetric. As a result of 
the sharp cut-off in the limb darkening, the kernel com- 
puted using the NextGen intensities in narrower than 
the black body kernel. However, for this particular situa- 
tion, the NextGen kernel has roughly the same width as 
the analytic kernel. Finally, at phase 180°, when the gi- 
ant is behind its unseen companion, all three kernels have 
more or less the same full width at half maximum. How- 
ever, the kernel computed using the black body intensities 
is flat-topped (Shahbaz ( 1998] ) also noted absorption line 
profiles with flat bottoms near this phase (his Fig. 1)), and 
the kernel computed with the NextGen intensities has a 
local minimum at zero velocity. The reason the kernels 
have flat tops or even central depressions near phase 180° 
is quite straightforward. The L\ point is in full view near 
this phase, and for high inclinations there are no points 
on the apparent stellar disk with fi = 1. For the exam- 



ple shown in Fig. ^, p, < 0.924 everywhere on the star at 
phase 180°. Hence the central parts of the apparent stel- 
lar disk are fainter (owing to limb darkening) than they 
would have been in the spherical case and there is less con- 
tribution to the kernel near zero velocity. If the intensities 
are computed using the NextGen table and the gravity 
is sufficiently low, then the central part of the star can be 
even darker. In some cases such as the one shown in Fig. 
^, the L\ point will be so dark that the broadening kernel 
will actually have a central depression near zero velocity. 



Given the variety of distortions in the broadening ker- 
nel as shown in Fig. |], it would be prudent to com- 
pute phase-specific broadening kernels when extracting 
V TO t sin i measurements from spectra (see also Shahbaz 



1998 ). As a specific example, we consider the eclipsing bi- 
nary UU Cancri, which consists of a K giant in a 96.7 day 
orbit about an essentially invisible (in the optical) com- 
panion (e.g. Eaton et al. |l99l| ). There is good evidence 
for an accretion disk around the unseen companion (Zola 
et al. 



1994) 



so we will assume that the K giant fills its 
Roche lobe exactly and that it is in synchronous rotation. 
Eaton et al. (1991) obtained a series of high resolution 
spectra and noted a change in the Doppler broadening of 
certain metallic absorption lines as a function of phase. 
This behaviour is consistent with expectations (Shahbaz 
1998). Using a linear limb darkening law with x = 0.85, 
Eaton et al. (1991) determined a rotational velocity of 
Vi-ot sin i = 25± 1 km s _1 for a spectrum near a quadrature 
phase and derived a mass ratio of Q = M comp /AI g i ant rs 
1.2 (using ELC, we derive a mass ratio Q = 1.235 from 
V ro t sini = 25 ± 1, where we have adopted a mass func- 
tion forthe unseen companion of f(M) — 0.56 M© (Pop- 
per 1977 )). We plot in Fig. |To| the analytic broadening 
kernel for x = 0.85 and a mean rotational velocity of 
2irR c ff/P — 25 km s , where R e s is the sphere-equivalent 
Roche lobe radius. We also show a broadening kernel com- 
puted using NextGen intensities, assuming Q = 1.65 and 
2ttR c s/P = 20.7 km s -1 . We have adopted an inclination 
of i = 89.6° and an effective temperature of 3900 K for 
the K giant (Zola et al. 1994). The two broadening ker- 
nels have nearly the same full width at half maximum. 
However, the implied component masses for the two cases 
are quite different. Assuming i « 90°, Q — 1.235 (the 
value from Eaton et al 

1.25 



M, 



1991) implies A/ g i ant 

logffgiant 



1.49 M Q , 
whereas our 



0.88 M Q , M, 



comp = 1.83 M Q , and 
value of Q — 1.65 gives M g i ant 
1.44 M©, and logg giant = 1.19. Zola et al. fll994| ) find a 
photometric mass ratio from their light curve solutions of 
q = M gian t/M comp — 0.564 ± 0.006, which in our notation 
is Q = 1.77 ±0.01. The case of UU Cnc is perhaps an ex- 
treme example, but it serves to illustrate the importance 
of using broadening kernels which account for deviations 
from linear limb darkening and spherical geometry. 
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Fig. 11. The curves show the integrated /-band specific 
intensity for various specific angles and effective tempera- 
tures as a function of the gravity. The intensity is generally 
a smooth function for a fixed temperature and emergent 
angle. Thus a linear interpolation in logg is reasonably 
accurate. 

3.3. Accuracy of the integration and interpolation 

Given the relatively large differences in the predicted light 
curves and broadening kernels that we predict, it is worth- 
while to discuss the numerical accuracy of the output light 
curves. Since the computation of the NextGen models 
is somewhat time consuming, we cannot tabulate mod- 
els for all effective temperatures and gravities. We have 
done extensive testing on the integration and interpola- 
tion procedures, and we believe the present table (steps 
of 0.5 dex in logg and 200 K in T e g) is a reasonable com- 
promise between CPU time and sampling accuracy. There 
are various ways to see how "smooth" the filter-integrated 
intensities are as a function of the effective temperature 
and gravity. Figs. |ll| and |l2| show two such representa- 
tions. In Fig. [ll] we show as a function of the gravity the 
integrated /-band intensity for four effective temperatures 
(4000, 4400, 5000, and 5600 K) at six different emergent 
angles (fj, = 0.05, 0.10, 0.15, 0.20, 0.25, and 0.30). In gen- 
eral, these curves are smooth, and a linear interpolation of 
the intensity in the logg direction gives reasonable accu- 
racy. Fig. |l2| shows a similar plot, but with the intensities 
as a function of the effective temperature. As before, the 
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Fig. 12. Similar to Fig. |lT], except the intensity is now 
displayed as a function of the effective temperature. Note 
the y-axis scale here is logarithmic. Within each panel, 
the specific angles are from bottom to top fi = 0.05, 0.10, 
0.15, 0.20, 0.25, and 0.30. In most cases the fi = 0.25 and 
the /i = 0.30 curves are nearly the same. As in Fig. O, 
the curves are quite smooth, and a linear interpolation of 
the intensity in T e g is reasonably accurate. 

curves are quite smooth and a linear interpolation in T e g- 
gives reasonable results. There is some irregular behaviour 
in the curves for fj, = 0.05 and fi — 0.10. In practice, how- 
ever, the contribution to the overall integrated light curve 
from angles less than 0.1 is small since the intensities are 
weighted by the value of fi. 

Another way to test the accuracy of the interpolation 
is to leave specific models out of the table, compute a light 
curve, and compare it to the light curve computed using 
the full intensity table. Fig. [l3] shows results of this ex- 
ercise for the star described in the caption of Fig. ||. We 
computed "regular" light curves where we used all mod- 
els in the intensity table and "cut" light curves where the 
T e g = 4000, logg = 1.0 model was excluded from the ta- 
ble. Thus, at T e ff = 4000, the local step size in logg was 1 
dex, rather than 0.5 dex. To compare the light curves, we 
define the "difference curve" as D = 100(/ r cg — ^cut)/Jreg ) 
where / rcg is the integrated light computed using the full 
intensity table and I cu t is the integrated light computed 
using the intensity table with the T e g = 4000, logg = 1.0 



12 



Orosz & Hauschildt: The use of NextGen model atmospheres in a light curve synthesis code 




Fig. 13. "Difference" plots in the light curves for the star 
described in the caption to Fig. [|. The top three sets of 
curves are difference curves for V, K, and /. For each filter, 
there are two curves corresponding to a "standard" grid 
of surface elements (N a = 40, Np = 14) and a "fine" grid 
(N a = 120, Np = 40). The curves for the fine grids are 
smoother. The curve at the bottom compares the regular 
V models computed using the two different grids. 



model removed. We see that the systematic difference be- 
tween the two curves is usually about 0.3 per cent and is 
at most 0.75 per cent in V. The corresponding values for 
/ and K are even smaller. We also computed light curves 
using two different grids of surface elements: the "stan- 
dard" grid with N a = 40 and Np = 14 and the "fine" grid 
with N a = 120 and Np = 40. Using the fine grid has little 
effect on the difference curves: the curves are smoother 
but the overall shapes are the same. The lowermost curve 
in Fig. |l^ compares the regular V light curve computed 
using the standard and fine grids. The two curves are the 
same to within 0.05 per cent (i.e. less than about 0.5 mil- 
limag). We conclude that for most situations our regular 
intensity table and the standard grid are adequate. 

There is nothing in our integration and interpolation 
techniques that limits us to « 1 per cent accuracy. Since 
our interpolation scheme does not have a fixed step size in 
the temperature or gravity we can easily add a few models 
with the appropriate temperatures and gravities for spe- 
cial cases where the data demand the highest accuracy. 
Other light curve codes also suffer systematic errors on 
the few per cent level. For example, most versions of the 
W-D code use a single limb darkening law for the entire 
star. Since the temperature and gravity are not constant 
over the surface of a tidally distorted star, small system- 
atic errors can be present when a limb darkening law for 
a single temperature and gravity are used. Van Hamme & 
Wilson (1994) have discussed this point in more detail. 



4. Discussion and summary 

In this paper we have presented a way to include specific 
intensities from detailed model atmosphere computations 
into a light curve synthesis code. We have shown that us- 
ing the model atmosphere intensities directly is almost re- 
quired for cool giants since the limb darkening behaviour 
for a cool low-gravity atmosphere is radically different 
than simple one- or two-parameter limb darkening laws 
(Fig. Q). This departure from the near-linear limb dark- 
ening law is a consequence of the spherical geometry used 
in the computation of the NextGen models. Other work- 
ers have computed spherical model atmospheres for cool 



giants before (e.g. Scholz & Tsuji 1984; Scholz & Takeda 
1992j ; J0rgensen et al. |1992j ), and the 



1987; Plez et al. 



strongly nonlinear limb darkening beh aviour has previ- 
ously been noted by Scholz & Takeda ( |l987[) , and more 
recently by Hofmann & Scholz (1998). We have shown that 
this strongly nonlinear limb darkening has a large effect 
on the light curve amplitude. In general, the light curves 
computed using the NextGen intensities have larger am- 
plitudes than the light curves computed using the same 
geometry but with black body intensities and a one- or 
two-parameter limb darkening law. Also, the strongly non- 
linear limb darkening has an effect on rotational broaden- 
ing kernels in that for spherical stars the kernel computed 
using the NextGen intensities will be narrower than a 
numerical or an analytic kernel that has a linear limb 
darkening law. If the star fills a substantial fraction of 
its Roche lobe, then the broadening kernel will also be 
different than the analytic one owing to the nonspherical 
shape of the star. 

Thus we basically have two different types of models 
with which to fit close binary light curves: our ELC code 
with the NextGen intensities, and codes that use either 
black body intensities or plane-parallel model atmosphere 
intensities and near-linear limb darkening (e.g. W-D, ELC 
in black body mode, ELC with the Kurucz table, etc.). It 
should be possible to test which class of models provides 
the better description of the observational data available 
for binaries with cool giants. In this regard, we have shown 
that observations in several bandpasses are useful in dis- 
criminating between the models. For example, Figs. [5] and 
[?] show how the difference in the depths of the minima 
between the NextGen and Kurucz light curves depends 
on the bandpass, while Fig. || shows that black body light 
curves provide a poor simultaneous fit to NextGen V 
and J light curves of a binary like T CrB. There is no 
shortage of potentially good binary stars to study. For ex- 
ample, T CrB, as we have already noted, has V, I, J, and 
K light curves available (although the light curves are per- 
haps somewhat noisy) . RZ Ophiuchi is an eclipsing binary 
consisting of a hot «B star and a K supergiant. The or- 
bital period is 262 days. Reasonably g ood l ight cu rves i n 
the bluer filters exist (e.g. Knee et al. 1986 ; Olson 1993 ), 
and additional observations in J, H , and K would be ex- 
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tremely valuable. Kruszewski & Semeniuk (1999) present 
a catalog of poorly studied eclipsing binaries with good 
parallaxes measured by the Hipparcos mission. Many of 
these binaries have long periods (more than 10 days) and 
contain evolved cool components. Needless to say, we en- 
courage observers to systematically study these binaries 
in multiple colours. With a binary of known distance, one 
has the added challenge of obtaining the correct integrated 
flux. 

Our interest should not be confined to binaries with a 
supergiant component. For example, there appears to be 
a small class of "long period" Algol type binaries (peri- 
ods between about 10 and 20 days) where the mass los- 
ing cool star appears to be slightly underfilling its Roche 
lobe. S ome examples are WW An dromed ae (Olson & Etzel 
1993a|) , S Can cri (O lson & Etzel |l993bj ), and DN Orionis 
(Etzel & Olson 1995 ). The cool stars in these three systems 
have surface gravities near logg ~ 2.5, so the sphericity 
effects are not as pronounced as they are in systems like T 
CrB. On the other hand, all three of these binaries have ex- 
cellent five-colour (7 {Kxon)ybvu) photometry, and S Cnc 
and DN Ori are totally eclipsing and are double-lined. Ol- 
son and Etzel report that in each case, the cool star seems 
to underfill its Roche lobe by about 10 per cent. It would 
be worthwhile to re-evaluate the light curve solutions for 
these three stars using ELC and the NextGen intensities 
to see if the sphericity effects can account for the apparent 
slight underfilling of the Roche lobes by the cool stars. 
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Appendix A: Brief outline of the ELC code 

A.l. Introduction and history 

Yorham Avni was interested in, among other things, mass 
determinations of the compact objects in high mass X- 
ray binaries (e.g. Cyg X-l, Cen X-3, etc.). During the 
course of his research he wrote a FORTRAN code to com- 
pute the ellipsoidal light curve of a single (usually) Roche- 
lobe filling star (Avni & Bahcall |1975| ; Avni |197Sfi . Avni 
passed on this code to Jeff McClintock and Ron Rcmillard 
shortly before his death in March of 1988 (see McClintock 



& Remillard 1990). The code was passed on to JAO from 
Ron Remillard sometime in 1994 when JAO was a gradu- 
ate student at Yale. During the spring of 1995 it became 
clear that the black hole binary GRO J1655-40 was an 
eclipsing system and that the Avni code in its original 
form was not adequate. The code was substantially modi- 
fied by JAO in the summer of 1995 to include light from an 
accretion disk and to account for eclipses (Orosz & Bailyn 



1997) 



The code described in Orosz & Bailyn (1997) is some- 
what unwieldy and it was becoming more and more dif- 
ficult to read and modify. Therefore a more general and 
more modular code was developed by JAO. Although we 
have used much of Avni's notation, and we have followed 
his basic method of setting up the Roche geometry and 
integrating the observable flux, most of the code is new. 
Owing to space limitations, we cannot give below each 
equation used. Rather, we give below some details of the 
parts of the code which have been substantially revised 
with respect to the earlier versions (there are numerous 
other papers and texts whic h go into var ying amounts 
of detai l, Wil son fc Dev inney |l97l|; W ilson |l979|; A vni & 
Bahcall |1975|: Avni h 9781: Linnell|l984|: Wilson Il990|: Orosz 



& Bailyn pL997| ) 



A. 2. The Potential 

We assume surface of the star is an equipotential surface of 
the following potential, which includes the gravitational, 



centrifugal, and Coriolis forces (Avni & Bahcall 1975) 



GM\ 
D 



1 m Q n 
1 Qx 

ri r 2 



1 + Q ( wi 



^orb 



(X 2 



(A.l) 



where M\ is the mass of the star under consideration, D 
is separation between two stars, Q = MijM\ is the mass 
ratio, uji is the star's rotational angular velocity, (jj or b is 
the orbital angular velocity, r\ and r 2 are the distances 
to the stellar centres in units of D, and x and y are nor- 
malised coordinates centred at star 1. For a binary with 
a given mass ratio Q, rotational angular velocity uj\ and 
orbital angular velocity a-> or b, there is a critical value of 
the potential Vtcrit where the star exactly fills its limiting 
lobe. Stars that are smaller then their limiting lobes will 
have * > * cri f 

To fully define the surface of the star, the user specifies 
Q, = wi/w orb , and the "filling factor" / = x po i B t/x L i, 
where a; p0 mt is the x-coordinate of the "nose" of the star 
and xli is the x-coordinate of the L\ point. In our case, 
the filling factor / is exactly 1 for Roche lobe-filling stars 
and less than 1 for detached stars. Situations with / > 1 
(the contact binaries) are currently not allowed. When 
/ < 1, the program computes x point from x point = fx L i 
and then computes ^(zpoint, 0, 0), which is the adopted 
potential for the star. Once the surface of the star is de- 
fined, a grid of surface elements is made using a polar 
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coordinate system with N a latitude rows equally spaced 
in the angle 9 and 4 * Np longitude points per latitude 
row equally spaced in the angle (f>. (Earlier versions of the 
code used a cylindrical coordinate system with N a rings 
equally spaced along the line of centres, running from the 
L\ point to the back of the star and 4 * Np surface points 
per ring equally spaced in angle.) It is convenient to use 
an internal rectangular coordinate system centred on star 
1 where the x axis points to the centre of the other object 
and the z axis is parallel to the direction of uj\. The value 
of \& and its derivatives are computed for each element, 
and from these quantities the local gravities g(x, y, z) and 
the surface normal vectors follow. 

A. 3. Mean temperature vs. polar temperature 

The temperature of the secondary star was defined in the 
Avni code by its polar temperature T po i e . Given T po i e , the 
temperatures of the other surface elements followed from 
the well-known von Zeipel relation: 



T(x,y,z) 



pole 



g{x,y,z) 



Spolc 



(A.2) 



The exponent (3 is 0.25 for stars with a radiative atmo- 
sphere (von Zeipel 1924) and 0.08 for stars with convec- 
tive envelopes (Lucy 1967). Wilson (1979) has pointed out 
that the polar temperature and the mean temperature of 
a distorted star will be different. In most cases, of course, 
one measures the mean temperature via the spectral type 



or colour index. Therefore, following Wilson (1979), we 
now have as input the mean (or equivalently effective) 
temperature of the star, denoted by T c s- The effective 
temperature is computed from the bolometric luminosity 



(A.3) 



where S is the surface area. T po i is then given by (Wilson 
1979) 



f pole 



/surface 9norm(x, V, z)^dS(x, y, z) 



1/4 



(A.4) 



where g n0 rm(x,y, z) = g(x,y,z)/ 'g po i c and dS(x,y,z) is an 
element of surf ace ar ea (Eqs. (A8), (A9), and (A10) of 
Orosz & Bailyn |l997j ). 



A.4- Addition of a second star 

Adding a second star to the code is relatively simple. We 
"flip" the mass ratio (define Q' — 1/Q), and solve for the 
potential and its gradients using the same subroutines as 
for the first star. When integrating the observed flux, we 
add 180° to the phase and use the same subroutines as for 
the first star. One complication occurs when the second 



star is not in synchronous rotation (Wilson 197E). In this 



case, the x derivative of the potential that is used in the 
detailed reflection routine (see below) has a different form. 



Thus the subroutine that returns the potential gradients 
returns two sets of x derivatives and the appropriate one 
is used for the detailed reflection. 



A. 5. Detailed reflection scheme 

Stars in a close binary can heat each other, and this mu- 
tual heating leads to easily observed consequences. Wilson 
( 1990] ) divides the "reflection" theory into four main parts: 
the geometric aspect, the bolometric energy exchange, the 
intensity from an irradiated stellar atmosphere, and the ef- 
fect on the envelope structure. In this paper Wilson gives 
a complete description of his "detailed reflection" scheme 
which treats the first two parts of the theory essentially 
exactly. We have fully implemented Wilson's treatment 
of the reflection effect, which is a big improvement over 



the scheme used by Orosz & Bailyn (1997). There are two 
points about this scheme that are noteworthy: 

(i) Wilson's scheme makes use of bolometric limb dark- 
ening approximations. We have seen that the limb darken- 
ing in cool giants is not well parameterised by the common 
limb darkening laws. Fortunately, in most practical situ- 
ations, a relatively hot star with a high surface gravity 
(log g > 4) irradiates a much cooler star. The irradiation 
of the hot star by the cool star can be neglected (in which 
case we don't care about the details of the cool star limb 
darkening), and the well-known limb darkening laws apply 
nicely to the hot, high-gravity star. 

(ii) Wilson introduced the use of the i?-function in the 
reflection scheme to allow for multiple reflection. For each 
element, the R- function is defined by 



Ri(x,y,z) = 1 
R 2 (x,y,z) 



F 1 



1 



Fi{x,y,z) 
F 2 (x,y,z) 



(A.5) 
(A.6) 



where F^/F^x, y, z) is the ratio of the total irradiating 
flux from star 2 seen at the local surface element on star 
1 and F[/F2(x, y, z) is the reverse. The new temperatures 
of the irradiated surface elements arc then 



Tr w (x,y lZ ) = Tr(x,y 1 z)R 1 1 /4 (x,y 1 z) 
TZ™(x,y,z) = T° ld (x, y ,z)R 1 2 /4 (x,y,z). 



(A.7) 
(A.8) 



It is usually assumed that the specific intensity of an ir- 
radiated surface element is the same as the intensity of 
an unheated surface element with the same temperature. 
At some point the irradiation become intense enough that 



this assumption must break down. Alencar & Vaz (1999) 



have computed some irradiated atmosphere models and 
presented limb darkening coefficients for use in light curve 
synthesis codes. The use of these coefficients is perhaps 
somewhat limited as the widely available versions of the 
popular W-D code have only a single limb darkening law 
for the entire star (which of course includes the unheated 
back hemisphere). One of us (PHH) is in the process of 
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computing irradiated atmospheres with the PHOENIX 
code. It should be relatively simple to include the irradi- 
ated atmospheres into the ELC code, at least for specific 
binaries. There exist many eclipsing close binaries where a 
hot subdwarf O/B star irradiates its G-M main sequence 
companion (e.g. Hilditch et al. 1996). Such systems can 



provide strong tests of irradiated atmospheres. 



A. 6. Flared accretion disk 



The accretion disk used in the Orosz & Bailyn (1997) code 
was a flattened cylinder. To make the disk perhaps more 
realistic, we have modified the disk so that its thickness 
as a function of the radius is proportional to the radius. 
The disk, if present, is always around star 2. Star 2 does 
not necessarily have to be present (as in an X-ray binary). 
The disk is described by five basic parameters: r outC r, the 
radius of the outer edge of the disk in terms of the effective 
Roche lobe radius of star 2; ri nncr , the inner radius of 
the disk in the same units as the filling factor of star 2 
(/i); Aim, the opening angle of the disk rim above the 
plane; Tdisk, the temp erature of the inner disk (in the 
Orosz & Bailyn ( 1997 ) code the temperature of the outer 
rim was specified); and £, the power-law exponent on the 
temperature profile of the disk: 



T(r) = T disk (r/r lnner )^ 



(A.9) 



For a steady-state accretion disk, £ = —3/4 (Pringlc 198l|) . 
For a disk heated by a central source, the exponent £ can 
take on a range of value s (—3 /4 ^ £ ^ — 1/2 (e.g. Fried- 
jung HH Vrtilek et al. |l99Cj ; Bell |1999| )). 

Since the surface of the disk is flared, each element on 
the face of the disk will have different "projection factors" 
r to the line of sight. Here T is equivalent to the angle 
/i discussed above. We define a polar coordinate system 
centred at the centre of the disk. The angle 9 is measured 
from the a;-axis in the direction of the positive y-axis. For 
a given the orbital phase cf>, T(r, 9) is given by 



F(r, 9) = — cos cj> sin i sin j3 cos 9 
— sin <j) sin i sin [3 sin 9 
+ cos i cos f5. 



(A.10) 



If T(r, 9) < 0, the point is not visible. 

There are several geometrical details which we must 
account for when the disk is flared and/or when star 2 
has a relatively large radius. First, for cases when the in- 
clination i is within (3 degrees of 90, parts of the disk face 
that have T(r, 9) > will be below the rim as seen by a 
distant observer. To account for these hidden points, we 
define the "horizon" of the top rim of the disk. In this 
discussion the horizon of an object is the outline of the 
object in sky coordinates. A point on the disk face is visi- 
ble if its sky coordinates are inside the top horizon. Star 2 
can block parts of the disk if its radius is relatively large. 
Since we currently require the inner radius of the disk to 



be equal to the radius of star 2 (if present) , some or all of 
the "bottom" part of star 2 may be hidden by the disk. 
Again, the shadowing of the disk by star 2 and the lower 
part of star 2 by the disk is easily accounted for by defining 
the suitable horizons. Finally, a disk with a large radius 
can inhibit the mutual irradiation of the two stars since 
the "top" of star 1 cannot "see" the "bottom" of star 2, 
and vice- versa. If a disk is present then inside the detailed 
reflection subroutine each line of sight between points on 
star 1 and star 2 is checked to see if it passes through the 
disk. 

Currently we assume that each surface element on the 
disk has a specific intensity that is the same as a nor- 
mal stellar atmosphere. Following Pringle et al. (1986) 



we use the model with the largest gravity for each effec- 
tive temperature. Of course, much more detailed model 
atmospheres specific to accretion disks are available (for 
example the grid of models for accretion disks in cata- 
clysmic variables presented by Wade & Hubeny ( |1998| )), 
and in principle a separate intensity table for accretion 
disks can almos t triv ially be added to ELC. Indeed, Lin- 
nell & Hubeny ( 1996 ) create light curves for binaries with 
disks by first computing detailed spectra for the disk using 
the Hubeny codes. However, the Hubeny code TLUSDISK 
is best suited for atmospheres hotter than about 10,000 K, 
and as such cannot really be used to model the cool outer 
parts of many disks. As it is now, our treatment of the ac- 
cretion disk is perhaps the most appropriate for systems 
where the disk is optically faint and its main effect on 
the light curves is geometrical (i.e. it eclipses the bright 
mass- losing star). Examples of such systems ar e W C rucis 
(Zola [1996D , GRO J1655-40 ( Orosz fc Bailyn |1997| ), and 
BG Gcminorum (Benson et al. 200C). 



A. 7. Third light 

In many Algol type binaries there is good evidence for a 
fainter third star that is gravitationally bound. For ex- 
ample, the O — C residuals of the eclipse timings of SW 
Lyncis are periodic and can be explained by the presence 
of a third body in a 5.8 year eccentric orbit about the 
inner binary (Ogloza et al. 1998). In such triple systems, 



the "third light" dilutes the observed amplitudes of the 
light curves from the binary, provided of course that the 
third star is sufficiently bright. We have a trivial way to 
self-consistently add third light to light curves in different 
bandpasses. In ELC, the user specifies three parameters 
for the third light: the temperature of the third star, its 
gravity, and its surface area relative to the surface area of 
star 1. The code interpolates the filter- integrated intensi- 
ties for the third star from the table, scales appropriately 
based on the surface area ratio, and adds the light to each 
light curve. 
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A. 8. Accuracy of light curves through eclipse 

The integration of the observed flux from a single star is 
straightforward and is sufficiently accurate for a reason- 
ably small number of surface elements. However, quan- 
tisation errors can become noticeable in the light curve 
of a star going through an eclipse. In many cases (e.g. 
grazing eclipses, no reflection effect) the number of sur- 
face elements can be modestly increased so that a smooth 
light curve can be obtained without a large increase in 
the CPU time required to compute the model. In other 
cases (e.g. deep eclipses and several iterations of detailed 
reflection) the number of surface elements needed to get 
smooth light curves becomes so large that the required 
CPU time becomes excessive. Thus ELC has two features 
which prove to be quite effective in greatly reducing the 
numerical noise of light curves though eclipse. 



A. 8.1. Improved horizon definition for the eclipsing star 

The horizon of the eclipsing star (the one "in front" ) is 
defined to be a collection of points on the star which have 
H = 0. In previous versions of the code, the program 
would step though the surface grid in the "a" direction 
and record which surface elements were last visible (i.e. 
the last point with fj, > 0). The resulting collection of 
sky coordinates of these surface elements would then de- 
fine the horizon of the star. There is a systematic error 
introduced when the number of surface elements is small 
since the numerical horizon of the star will be slightly 
smaller than the actual horizon. In the current version of 
the code, the program steps along each latitude row and 
records the ^-coordinate of the last visible point V ; S and 
the (^-coordinate of the first point hidden below the hori- 
zon 4>hid- A simple bisection procedure is used to find the 
(^-coordinate (with a given latitude #;) where \x — 0. Fif- 
teen iterations of this bisection procedure are enough to 
find ^hor to better than 10~ 5 radians when Np — 6 (i.e. 
24 longitude points). The corresponding angles [i are all 
^5 5 x 10~ 6 . The x and y sky coordinates of the point with 
the surface coordinates #i,</>hid are then determined. A 
similar procedure is done where the program steps through 
the 9 angle for each longitude <f>i and #hor is found from 
v is and #hid using bisection. After the list of points on the 
star with [i — is generated, the x and y rectangular co- 
ordinates of each point on the sky are converted to a i?,0 
polar coordinate system and sorted in the polar angle 0. 
The sorted array forms a convex polygon on the sky. If 
the number of surface elements on the star is relatively 
small (N a ^ 30 and Np ^ 6), then the actual horizon of 
the star can have some curvature between adjacent points 
on the polygon. Since the radius R as a function of is 
always a very smooth function, we use spline interpola- 
tion to resample R for every 1° in 0. The new resampled 
polygon with 360 points always has enough points so that 
the horizon of the star essentially has no curvature be- 



tween adjacent points. We have done numerous tests and 
found that the horizons derived using a small number of 
surface elements with the new routine are always the same 
as the horizons found by the old routine with a very large 
number of surface elements (N a ss 400 and Np w 100). 

A. 8. 2. Fractional surface elements on the eclipsed star 

To compute the observed flux in a given bandpass from a 
star at a given phase, we numerically evaluate the integral 
in Eq. || using all of the surface elements with /i > 0: 

hlLT = ^2 ^2 lFIl,T(fJ'i,j)(-li,jrl j A0A(t>/ COS Pi j (A. 11) 

i=l j=l 

where AO and A<j) are the angular spacings of the ele- 
ments in latitude and longitude, respectively, and where 
(3 is the angle between the surface normal and the ra- 
dius from the centre of the star. In other words, we per- 
form a simple numerical quadrature along each latitude 
row where the points are equally spaced in the angle <fi. If 
there is an eclipsing body in front of the star whose flux 
is being evaluated, then each point on the star in back is 
projected onto the sky and a simple routine is used to see 
if this point is inside the polygon representing the horizon 
of the body in front. If the point in question is eclipsed, 
its flux contribution is simply left out of the summation. 

In general, the horizon of the star in front will not pass 
exactly between two adjacent points on a given latitude 
row on the star in back. As a result there will be "visi- 
ble" points that are actually centred in partially eclipsed 
surface elements and "eclipsed" points that are centred in 
partially visible surface elements. If the number of surface 
elements on the star in back is large enough, the contri- 
bution from the fractionally eclipsed surface elements will 
tend to cancel out. However, it is much more computa- 
tionally efficient to use a smaller number of surface ele- 
ments and make a correction for the fractionally eclipsed 
pixels. At each latitude row 6i on the star in back, the pro- 
gram determines the (^-coordinate of the last point visible 
before the horizon of the eclipsing body v is and the (j>- 
coordinate of the first point hidden behind the horizon 
of the eclipsing body 0hid- Another bisection procedure is 
used to determine the 0-coordinate on the star in back 
where the horizon of the eclipsing body intersects the 0i 
latitude row </f>hor- If l^hor — 4>vi&\ < A<p/2, then the last 
point visible before the horizon is centred in a partially 
eclipsed surface element and a negative correction is added 
to the flux summation: 

2(0hor - 0vis) — A(j>] 



TP — 
± r.orr — 



A(j) 



where F v - ls is flux from the last visible point, i.e. 



F V is — /filt (Mi,i)M« 
Likewise, if |^>hor — 



jAOAcf)/ cos/3, 



(A.12) 



(A.13) 



0hid| < A(f>/2, then the first point 
hidden behind the horizon is centred in a partially visible 
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surface element and a positive correction is added to the 
flux summation: 



Fr, 



A(/) - 2(</>hor - </>hid) 



(A.14) 



where Fhid is the flux the eclipsed point. If there is an 
annular eclipse, then this procedure for finding </>h r may 
have to be done twice for a given latitude row. We have 
tested this simple interpolation procedure quite exten- 
sively and have found it to be quite effective. Smooth light 
curves through eclipse can be obtained for grid sizes as 
small as N a = 14 and Np — 6. 

A. 9. Comparison with Wilson and Devinney 

We have done extensive testing of the ELC code in its 
black body mode against the W-D code. Unfortunately, 
there is some confusion over the notation of some of the 
input parameters between the two codes. In particular, f2 
in ELC is the ratio of the star's angular velocity to the 
orbital velocity, whereas in W-D refers to the potential. 
The W-D il-potentials essentially define the shapes of the 
stars. Perhaps to conserve parity, / in ELC is the "filling 
factor" which has the same function as the Sl-potential 
in W-D (i.e. it defines the shape of the star), whereas in 
W-D, F is the ratio of the star's angular velocity to the 
orbital velocity (ELC's f2). The phase convention between 
the two codes is different. In ELC, star 1 is in front of star 
2 at phase 0.0, whereas in W-D star 1 is behind star 2 at 
phase 0.0. 

To facilitate comparisons between ELC and W-D, we 
do two things. First, since the internal form of the po- 
tential is the same for both codes, ELC prints out the 
program values of the potentials which then are the input 
ft potentials for W-D. Second, we add a phase shift of 0.5 
to the W-D light curves (the 'pshift' input parameter). 
We computed various model binary light curves using the 
two codes and compared them by normalising the light 
curves at phase 0.25 (quadrature). We of course used ex- 
actly the same effective wavelengths, limb darkening laws 
and coefficients, and reflection schemes for both codes. In 
most cases, the light curves agreed to better than 0.1 per 
cent (better than 1 millimag). If we compared the light 
curves by adjusting the normalisation of one of the curves 
to match the other, then the largest deviations became 
even smaller. 

ELC also computes radial velocity curves (Wilson & 
Sofia 1976). The velocity curves from ELC agree with the 



W-D velocity curves to better than 0.1 per cent. 

Finally, we compared computed geometric quantities 
between the two codes (i.e. the "polar" radius, "point" ra- 
dius, etc.). The agreement was essentially exact. We also 
computed sphere-equivalent Roch e lobe radii and com- 
pared our results with Eggleton's ( 1983 ) results and like- 
wise found nearly exact agreement. 
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